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Summary. We describe a family of quadrilateral meshes based on diamonds, 
I ^ rhombi with 60° and 120° angles, and kites with 60°, 90°, and 120° angles, that 

can be adapted to a local size function by local subdivision operations. The vertices 
of our meshes form the centers of the circles in a pair of dual circle packings in which 
each tangency between two circles is crossed orthogonally by a tangency between 
two dual circles. 
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1 Introduction 



A famous and deep theorem of Koebe, Andreev, and Thurston [Tj[2{ |23||40] 
I asserts that the vertices of every planar graph may be represented by a circle 

packing, a system of circles with disjoint interiors, such that two vertices are 
CN| adjacent in the graph if and only if the corresponding two circles are tangent. 

This representation is not unique without additional constraints (for instance, 
a 4-cycle has infinitely many distinct representations as a set of four tangent 
circles) but it can be made unique, up to Mobius transformations, in one of 
two different ways: 

O^l • Let G be constrained to be a maximal planar graph; that is, every face 

of G, including the outer face, must be a triangle. Then its representation 
as a circle packing exists and is unique up to Mobius transformations 
(Thurston, Corollary 13.6.2). We call this a maximal circle packing. An 
example is shown in Figure [I] left. 
• Alternatively, let G be a 3-vertex-connected planar graph. It has a unique 
planar embedding; let G" be the dual graph of this embedding. Then it is 
possible to represent both G and G' by simultaneous circle packings with 
the property that, for every edge e of G and its corresponding dual edge 
e', the two circles representing the endpoints of e have the same point of 
tangency as the two circles representing the endpoints of e' and, moreover, 
the circles for e cross the circles for e' at right angles at this point. Again, 
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Fig. 1. A maximal circle packing (left) and an orthogonal circle packing (right). 



this representation is unique up to Mobius transformations [10] , and we 
call it an orthogonal circle packing. An example is shown in Figure [T] right. 

Circle packings of these types have been applied in the field of graph draw- 
ing, to find drawings of planar graphs with right angle crossings [10], high 



angular resolution 13 29 , and small numbers of distinct slopes [21] . They 
have also been used for mesh partitioning 17 32p3 , for visualization of brain 



structures 20 , for analyzing the structure of soap bubbles ]16], for solving 
differential equations [19] , for constructing Riemann surfaces from combinato- 
rial data [9J, and for finding approximations to conformal mappings between 
different simply connected domains, which can be used as an important step 
in structured mesh generation [34 , 37) . 

The two constrained forms of circle packing guaranteed to exist by the 
circle packing theorem would seem to be also a natural fit for unstructured 
mesh generation: in a maximal circle packing, the graph of adjacencies be- 
tween tangent circles (with its vertices placed at the triangle centers) forms 
an unstructured triangle mesh, and in an orthogonal circle packing, the graph 
of adjacencies between orthogonal circles forms an unstructured quadrilateral 
mesh. Additionally, if the degree of a graph is bounded, then the circle pack- 
ings generated from it are naturally graded in size: adjacent circles have radii 
whose ratio is bounded, and the triangular or quadrilateral elements derived 
from the packing have bounded aspect ratio. However, despite their obvious 
appeal, these types of circle packing have not been used in mesh generation, 
because the geometry of a circle packing is difficult to control: circle packings 
are generated from combinatorial data (a graph) rather than from geometric 
data (the shape of a domain to be meshed) and in general, a small localized 
change to the graph from which the circle packing is generated can lead to 
large and non-localized changes to the packing. 

Circle packings are still used in many unstructured mesh generation al- 
gorithms, but they are not the types of packings described by the Koebe- 
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Andrccv-Thurston theorem. Rather, methods for unstructured mesh genera- 
tion based on circle packings have typically relied on methods that generate 
packings with less structure, by placing circles one at a time using a greedy 
algorithm [4j[7j[T5||24j[25j[27j|4l], physical simulation 36 , or decimation of 



quadtree-based ovcrsampling schemes 30 31 . Circle packing mesh generation 



algorithms have been used to find nonobtuse triangular meshes for polygonal 
domains [7|[l5j as well as bounded-aspect-ratio triangular meshes [25] and 
quadrilateral meshes in which all elements belong to certain special types of 
quadrilateral p] . The circle packings generated by these methods can be made 



to have radii controlled by a local size function 25 41 , and mesh generation 



techniques based on these methods can be applied in higher dimensions as 



well 24 28 30 31 . However, these circle packings are neither maximal nor 
orthogonal; rather, as well as having the three-sided gaps between circles that 
a maximal circle packing would have, they also include irregular gaps with 
four or more sides. Instead of generating a mesh directly from the intersection 
pattern of the given circles, these methods need additional vertices at points 
such as the circumcenters of the gaps, and they typically also need additional 
case analysis to handle the different possible shapes of their gaps. 

In this paper, we show for the first time that it is possible to construct or- 
thogonal circle packings, one of the two types of circle packing guaranteed to 
exist by the Koebe-Andreev-Thurston theorem, with geometric rather than 
graph-theoretic control of the position and size of the circles. Our circle pack- 
ings are based on a quadrilateral mesh in which every quadrilateral is either 
a diamond, a rhombus with 60° and 120° angles, or a kite with 60°, 90°, 
and 120° angles; the circles in the packing are centered at the vertices of this 
mesh and have their points of tangency at the points within each quadrilateral 
where its two diagonals cross. 

In some ways, our construction resembles more standard meshing tech- 
niques based on quadtrees (squares recursively subdivided into four smaller 
squares) [5]|T4p8p2] , or the recursive subdivision of triangles into four smaller 



triangles used by the Sloan digital sky survey 39 . Like quadtrees and recur- 



sive triangle meshes, any two quadrilateral meshes formed by our construction 
can be transformed into each other by a sequence of simple local operations, 
and these local operations give the set of meshes formed in this way the struc- 
ture of a distributive lattice, in which any two meshes have a coarsest common 
refinement and a finest common coarsening. Again, like quadtrees, these sub- 
division operations allow our mesh to be adapted to a size function specifying 
the element size at each point of a domain, and meshes adapted to a size 
function in this way have a number of elements that is within a constant 
factor of optimal for any bounded-aspect-ratio mesh for that size function. 
Like quadtrees, our meshes can as well use the same subdivision operations to 
adapt dynamically to varying size functions at different time steps of a finite 
element simulation. However, the squares of a balanced quadtree may share 
a boundary edge with up to eight smaller squares, and therefore require ad- 
ditional subdivision to produce a triangle or quadrilateral mesh; in contrast, 
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Fig. 2. Left: The rhombille tiling. Right: The short diagonals of the rhombille tiling 
form a hexagonal tiling. 



our method provides a mesh directly without the need for additional sub- 
division. The hexagon-based meshing algorithms of Sufincr, Grciner, Liang, 
and Zhang 26 38 also resemble the method here, but are based on scaling 
by factors of two at each level of subdivision whereas we use factors of v3- 
Other subdivision schemes related to the methods described here include hon- 

and \/3 refinement 22 43 ; however, these schemes 
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eycomb refinement 

are typically applied to every element of a structured or semistructured mesh 
rather than (as here) to selected elements of an unstructured mesh. 

We caution that our results should be viewed as primarily mathematical 
rather than being ready for immediate use in mesh generation practice. The 
meshes described here do not conform to domain boundaries with arbitrary 
slopes, and we have not implemented our algorithms. 



2 The rhombille tiling and its local subdivision 

Our construction begins with the rhombille tiling of Figure [2] (left), a tessel- 
lation of the plane by rhombi with 60° and 120° angles that resembles the 
axonometric projection of a pile of cubes and that can be formed by subdi- 
viding a regular tiling by hexagons into three rhombi per hexagon [11| . Each 
vertex of the rhombille tiling has degree (valency) either three or six: either 
six rhombi meet at their acute corners, or three rhombi meet at their obtuse 
corners. The short diagonals of the rhombi form another hexagonal tiling in 
which each hexagon surrounds a degree-six vertex of the rhombille tiling. 

Suppose that we wish to refine a rhombille tiling within a region R of 
the plane, forming a mesh with smaller elements, while leaving the tiling un- 
changed outside R. To do so, we may approximate R by a set of the hexagons 
formed by short diagonals, and perform the local replacement operation il- 
lustrated in Figure [3] within each hexagon. This operation replaces the six 
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Fig. 3. Replacement of six rhombille-tiling edges within a hexagon of short diagonals 
(left) by six rhombi with side length l/y/3 times the lengths of the sides in the 
original tiling. 




Fig. 4. The result of performing multiple local replacements in a rhombille tiling. 
Within the blue replaced region, we have another rhombille tiling, rotated from the 
original and with smaller tiles. 

edges that meet in the center of the hexagon with a network of 18 edges, 
shorter by a factor of l/-\/3 from the original edges. These new edges form the 
boundaries of six rhombi similar to the ones in the original rhombille tiling 
but rotated from them by 30° angles. Each subdivided quadrilateral crossing 
the boundary of the hexagon remains a quadrilateral after the replacement, so 
the result of the replacement is a valid quadrilateral mesh with six additional 
quadrilaterals for every replaced hexagon. 

Figure [4] shows the mesh resulting from the performance of multiple local 
replacements in a rhombille tiling. When one of the replaced hexagons abuts 
a hexagon that has not been replaced, the quadrilaterals formed across the 
shared boundary of the two hexagons are kite-shaped, with vertex angles 60° , 
90°, and 120°. However, when two or more replaced hexagons abut each other, 
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Fig. 5. Changes to the system of circular arcs within each quadrilateral caused by 
a single local replacement step. 

the quadrilaterals that lie across their shared boundary are rhombi congruent 
to the ones contained within each replaced hexagon. Within any region formed 
by multiple replaced hexagons, the smaller rhombi formed by this replacement 
process meet up with each other in the pattern of another rhombille tiling, 
rotated from the original tiling by a 30° angle and with tiles that are smaller 
by a 1/V3 factor. 

Once this replacement has been performed, the same replacement process 
can be performed within the finer rhombille tiling formed within the replaced 
region. The degree-six vertices that can be replaced within this finer tiling 
are either the same degree-six vertices that were replaced previously, or the 
points where three replaced hexagons meet. 

We call the meshes generated by any number of steps of this replacement 
process "diamond-kite meshes" . 

3 Orthogonal circle packing 

We now show that the diamond-kite meshes correspond to orthogonal circle 
packings. In each quadrilateral of a diamond-kite mesh, place arcs of four cir- 
cles, centered at the quadrilateral's four vertices and meeting at the center 
of the quadrilateral. Then, the circular arcs for the quadrilaterals meeting at 
a vertex will necessarily link up to form a single circle. In the unsubdivided 
rhombille tiling, this is true because the quadrilaterals sharing a vertex are all 
rotated images of each other, and it remains true in each of the local replace- 
ment steps by which the subdivided tiling is formed. As shown in Figure[5j the 
circular arcs surrounding each vertex of the replaced hexagon (shown as green 
in the figure) retain their previous radius, and the circular arcs surrounding 
the center vertex and each newly added vertex (shown as violet) meet up to 
form a circle that lies entirely within the replacement region. 

The circles formed in this way meet in tangent pairs at the points within 
each quadrilateral where the diagonals cross, and (because the rhombs and 
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Fig. 6. A diamond-kite mesh with quadrilaterals of many different scales and its 
circle packing. 

kites of the mesh are both orthodiagonal) the two pairs of tangent circles 
meeting at each crossing point are orthogonal to each other. Thus, the result 
is an orthogonal circle packing. Figure [6] shows a larger example. 

4 Lattice of refinements 

Given an initial rhombille tiling T, we may uniquely specify each of the local 
replacement steps used to form a diamond-kite mesh by a pair of parameters 
(p, s) , where p is the center point of the replaced hexagon and s is its side 
length. We may analyze cases to determine the conditions under which a 
replacement with parameters (p, s) is possible: 

• If s is the length of the sides of the original tiles of T, then replacement 
(p, s) may always be performed in every diamond-kite mesh formed from 
T in which it has not already been performed. 

• If s is smaller than the side length in T, and (p, sy/3) is the pair of pa- 
rameters for another replacement step, then replacement (p, s) may be 
performed if and only if replacement (p, sy/3) has already been performed. 
The reason for this is because the replacement (p, sy/3) is the only possible 
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way to incorporate into the tiling the six edges that will be replaced by 
(P, s). 

• In the remaining case, there are three points po, pi, and P2 equally spaced 
around p at distance sv3 from it, such that each point pi gives the pa- 
rameterization of a replacement step (pi, sy/3) and such that replacement 
(p, s) may be performed if and only if all three of (pi, sV3) have already 
been performed. Each of the three replacements (pt, sv3) is the only pos- 
sible way of incorporating into the tiling two of the six edges that will be 
replaced by (p, s) . 

We may encode this case analysis by an infinite directed acyclic graph Gt 
in which we have a vertex for each pair (p, s) that defines a valid replacement 
step, and in which we have an edge from (p, s) to (p, sy3) (when (p, sa/3) 
defines a valid replacement step) or from (p, s) to each of the three vertices 
(pi,sy/3) (in the cases where these three vertices are defined). In Gt, the 
replacement steps corresponding to all incoming neighbors of a vertex must 
be performed before the replacement step corresponding to the vertex itself 
may be performed. 

Alternatively and equivalently, we may use an infinite partially ordered 
set Pt that has an element for each pair (p, s), and in which two elements 
x and y are ordered x < y whenever there is a directed path from x to y in 
Gt- This formulation is more convenient for describing the set {x \ x < y} of 
all replacement steps that must be performed prior to performing operation 
y, either because they directly precede y (as in the case analysis above) or 
because they precede one of the predecessors of y. 

If T' is a diamond-kite mesh formed by refining the initial mesh T, then 
(by the analysis above) the set of replacement steps that were used to generate 
T' from T must be a finite lower set in Pt, that is, a set L of elements with 
the property that, for every element y £ L and for every element x with x < y, 
x is also in L. Conversely, every finite lower set L uniquely defines a diamond- 
kite mesh Tl as a refinement of T, because the replacement operations in 
L may be performed in any order that is consistent with the case analysis 
above; different orderings of the same operations will always lead to the same 
results. The orderings in which a given lower set L of replacement steps may 
be performed are exactly the linear extensions of the partial order induced by 
the elements of L in Pt, and a valid ordering for a given lower set L may be 
calculated by applying a topological sorting algorithm to the directed acyclic 
graph induced as a subgraph of Gt by the elements of L. 

Let Ti and Ti be any two diamond-kite meshes formed by refining T, and 
described by the respective finite lower sets Li and L 2 in Pt- Then the sets 
Li n L 2 and L\ U Li are also finite lower sets, describing respectively the finest 
mesh Ti A Ti from which both T\ and T2 can be formed by refinement, and 
the coarsest mesh T\ V Ti that can be formed by refining both T\ and T2. In 
this way, as with the lower sets of every partially ordered set [8] , the family 
of diamond-kite meshes can be given the structure of a distributive lattice. 
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Fig. 7. Cases for refine(p). From left to right: (a) p has degree six, and is ready for 
immediate replacement; (b) p has degree three, and does not meet the preconditions 
of refine; (c) p has degree five, and a replacement at the 60° kite vertex must be 
performed prior to a replacement at p; and (d) p has degree four, and two 60° kite 
vertices must be replaced prior to a replacement at p. 

We remark that the same ideas of constructing an infinite graph describ- 
ing the prerequisite relation between potential replacement steps, deriving an 
infinite partial order from the graph, and describing each possible mesh as a 
lower set of this partial order, can be applied equally well to describe the set 
of possible balanced quadtrees derived from an initial square. Thus, the set of 
balanced quadtrees can also be given the structure of a distributive lattice. 

5 Local replacement at mesh vertices 

When we adapt a mesh to a local size function, it will be convenient to specify 
the local replacement steps of the adaption process by vertices of a mesh 
rather than by the more abstract elements of an infinite partially ordered set 
described in the previous section. 

We will define a recursive subdivision algorithm refine (p) that performs 
a single local replacement step at a vertex p of a diamond-kite mesh, after 
performing all prerequisite replacements, as follows. We require, as a precon- 
dition for this algorithm, that p be a 60° vertex of at least one mesh element. 
Based on this requirement, the case analysis below describes which other re- 
placement steps are necessary before the one at p can be performed: 

• If p has degree six, then the most recent replacement step that affected 
the edges incident to p must either have replaced a hexagon with p at its 
center, or have been the third of three replacements of hexagons meeting 
at p. In this case, p is already surrounded by diamonds and/or kites having 
60° angles at p, as shown in Figure [7^ a) . It is possible to perform a local 
replacement step at p without performing any other replacements. 

• If o has degree three, then it must be the case that the most recent re- 
placement step at p replaced a hexagon for which p was interior but not 
central. Then p is a 120° vertex of three elements, either two diamonds 
and one larger kite (as in Figure ufb)) or three diamonds. In this case, 
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it is not possible for the precondition of the refine algorithm to be met, 
because there is no 60° angle at p. 

• If p has degree five, then it must be the case that the most recent replace- 
ment to affect the neighborhood of p was the replacement of a hexagon 
having p as a vertex, and additionally this must have been the second 
replacement of the three hexagons of that size meeting at p. Then the 
neighborhood of p consists of three elements with 60° angles (diamonds 
or kites within the two replaced hexagons) and two elements with 90° an- 
gles (necessarily kites in the third hexagon); see Figure [7^c). Let q be the 
shared 60° vertex of these two kites. We recursively call refine(g), after 
which we may perform a replacement step at p. 

• If p has degree four, then p was a vertex of the hexagon for the most recent 
replacement at p, and this was the first replacement of the three hexagons 
of that size meeting at p. In this case p has a neighorhood with one 60° 
angle (a diamond or kite in the replaced hexagon), two 90° angles (kites 
in the two unreplaced hexagons), and one 120° angle (a rhomb or kite 
overlapping the two unreplaced hexagons), as is depicted in Figure J7jd). 
We recursively call refine for each of the two 60° vertices of the kites with 
incident 90° angles, after which we may perform a replacement step at p. 

Thus, refine(p) need only list all kites that have 90° angles at p, recursively 
call itself on the 60° angles of these kites, and then perform a replacement step 
at p. The recursion necessarily terminates, because each recursive call leads 
to a replacement step on a larger hexagon. Each recursive call adds elements 
to the mesh for the replacement step it performs, so the total time for this 
recursive procedure is linear in the total change to the number of elements in 
the mesh. The case analysis above shows that this recursion performs exactly 
the replacement steps that are predecessors of (p, s) in the partial order P$ 
and that had not already been performed at the start of the recursion. 

6 Adaption to a local size function 

We define a local size function to be a function a that maps each point p 
of the plane (or of a subset of the plane to be meshed) to a positive real 
number <t(j>), specifying the largest allowable side length of a mesh element 
containing p^ We assume that access to a is via a subroutine oversized(Q) 
that takes as argument a quadrilateral Q and returns a Boolean value, true if 
Q contains a point p for which a(p) is less than the side length of Q and false 
otherwise. Our task is to find a mesh that is as coarse as possible subject to 
the constraint that oversized returns false for all mesh elements. 

It would be equivalent to within constant factors to specify the maximum allow- 
able area, perimeter, diameter, or circumradius of the element, but side length turns 
out to be more convenient for our purposes, because it leads to fewer ambiguities 
about which replacement steps are necessary. 
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To do so, from an initial mesh, we make a queue of unprocessed quadri- 
laterals in the mesh. We repeatedly find and remove a quadrilateral Q from 
this queue, process it, and add to the queue all new quadrilaterals formed 
during processing, until the queue becomes empty. To process a quadrilateral 
Q, we first try calling oversized(Q); if it returns false, processing is complete. 
Otherwise, if Q is a kite, we let p be its 60° vertex and call refme(p). If Q is 
a diamond, let pi and pi be its 60° vertices and let Qi an d Qi be the kites 
within Q that have the same maximum side length as Q and that have pi 
and P2 (respectively) as their 60° vertices. For each i in the set {1, 2} we call 
oversized(Qi) and, if it returns true, we call refine (pi). Note that it is not pos- 
sible for both oversized(Qi) and oversized (Q 2) to be false, because that would 
be inconsistent with the assumed true value of oversized(Q). Thus, regardless 
of the shape of Q, if oversized(Q) is true then processing Q will cause it to 
become subdivided. 

Whenever this adaption procedure calls refine (p), leading to a replace- 
ment step at vertex p, the same replacement step must be performed in all 
diamond-kite meshes that obey the size function <j, because otherwise the 
current quadrilateral Q or a kite within it would remain and have too large 
a side length. Therefore, the result of the adaption procedure is the coarsest 
diamond-kite mesh consistent with the size function. Each step either removes 
a quadrilateral from the queue without adding any others, or it takes time 
linear in the number of elements added, so the total time for this adaption 
procedure is linear in the size of the final mesh it produces. 



7 Size and length optimality 

Following Ruppert |35| , we may use local feature size to prove that the method 
for size adaption discussed in the previous section produces meshes that (com- 
pared to any other mesh with quadrilaterals or triangles of bounded aspect 
ratio obeying the constraints of the local size function) are within a constant 
factor of optimal with respect both to their number of elements and to their 
total edge length, matching known results for quadtree meshes [5 14 and for 
meshes formed by Delaunay refinement [35] . 

We assume that the size function cr(p) is defined in such a way as to lead 
to a finite mesh, and we define the local feature size to be a function cr(p) that 
maps a point p to the number 

<r(p) = inf {distance(p, q) + o-(q)}. 

The point q in the minimization ranges over the rest of the plane, but its 
minimum will necessarily occur within a disk of radius cr(p) centered at p, 
because all other points lead to larger values than the value <r(p) achieved at 
q = p. The following two observations are central to our analysis: 

• In all meshes with bounded-aspect-ratio elements, the size of the element 
containing a given point p is f2(<j(p)). To see this, let q be a point achieving 
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(or approximately achieving) the minimum value in the definition of <r(p). 
The element containing q must have size at most cr{p), and the sequence 
of elements crossed by the line segment from q to p cannot increase in size 
from that value by more than a constant factor before they reach p. 
• In the diamond-kite mesh defined from the size function a, the size of 
the element containing p is 0(a(p)). This follows from the fact that, if q 
is any point in the plane, the size of the smallest element of the partial 
order Px that is forced by the value of cr(q) to be included in the mesh 
and that corresponds to a local replacement for a hexagon containing p is 
0(distance(p, q) + <r{q). 

Given a bounded-aspect-ratio mesh M, let a(p) denote the area of the 
element containing p. Then, for all elements E of M, we have the identity 

Je ol{p) 

Therefore, the number of elements in the mesh can be counted by 

dp. 



M Ot(p) 

However, l/a(p) is lower-bounded by fl((a(p))~ 2 ) for all bounded-aspect- 
ratio meshes, and upper-bounded by 0((<j(p))~ 2 ) for diamond-kite meshes; 
therefore, the number of elements in the diamond-kite mesh determined by 
size function a is within a constant factor of optimal. 

Similarly, given a bounded-aspect-ratio mesh M, let 7r(p) denote the 
perimeter of the element containing p. Then, for all elements E of M, we 
have the identity 

tt( P ) 



L 



. dp = perimeter (p). 
e a{p) 

Therefore, the total perimeter of all the elements in the mesh is 

Jm a (p) P ' 

However, ir(p)/a(p) is lower-bounded by f2(l/a(p)) for all bounded-aspect- 
ratio meshes, and upper-bounded by 0(l/a(p)) for diamond-kite meshes; 
therefore, the total perimeter of the elements in the diamond-kite mesh de- 
termined by size function a is within a constant factor of optimal. 



8 Additional properties 

Every vertex of a diamond-kite mesh has degree at most six. Every two or- 
thogonal circles have radii differing by a factor of exactly and every two 
tangent circles have radii differing by a factor of at most three. 
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Fig. 8. Left: A 3-colored diamond-kite mesh. Right: A quadrilateral mesh that 
requires four colors in any face coloring. 

Graph colorings of meshes may be used to schedule batches of parallel 
updates to the values stored at mesh elements, in order to ensure that each 
two values that are updated in the same batch are independent from each 
other [3j[6] . As with all planar graphs in which every face has an even number 
of sides, the vertices of a diamond-kite mesh may be colored with two colors, 
but in this context it is more relevant to color the faces of the mesh so that no 
two faces that share an edge have the same color. This may be done with three 
colors by the following simple strategy: define an equivalence relation on the 
quadrilaterals of the mesh, according to which two quadrilaterals are equiva- 
lent when their diagonals are parallel, and assign one color to each equivalence 
class. There are only three equivalence classes: quadrilaterals in two different 
equivalence classes will have their diagonals rotated by 30° from each other, 
and after three such rotations we return to the starting equivalence class. No 
two adjacent quadrilaterals in a diamond-kite mesh may have parallel diago- 
nals, so adjacent quadrilaterals are always assigned distinct colors. Therefore, 
the result is a proper 3-coloring. In contrast, quadrilateral meshes other than 
the diamond-kite meshes may sometimes require four colors (Figure [8]). 

9 Conclusions 

We have defined the family of diamond-kite meshes based on a simple lo- 
cal replacement step starting with the rhombille tiling. In these meshes, the 
most acute angle is 60°, the most obtuse angle is 120°, and all elements have 
bounded aspect ratio. The element size can be controlled by a local size func- 
tion, and the number of elements and total edge length of the elements is 
within a constant factor of optimal for the given size function. Replacement 
operations may be performed adaptively to handle time-dependent size func- 
tions. Unlike adaptive quadtree meshes, this system provides a quadrilateral 
mesh directly without any need for additional subdivision. 
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The most novel feature of our new meshes is that their vertices form the 
centers of an orthogonal circle packing of the type guaranteed to exist by the 
Koebe-Andreev-Thurston circle packing theorem, showing for the first time 
that it is possible to incorporate this type of circle packing into a meshing 
algorithm. 

Much remains to be studied in this area. On the mathematical side, we still 
do not know whether it is possible to define an analogous local replacement 
scheme that would allow the generation of maximal circle packings (in which 
every gap between three circles has exactly three sides) with similar properties 
to those of the diamond-kite mesh, and in particular with the ability to control 
the size of the circles in one part of the packing without changing the geometry 
of the circles in distant parts of the packing. On the more practical side, it 
would be of interest to develop the diamond-kite method into a practical mesh 
generation system and to compare empirically the quality of meshes generated 
in this way with those from other comparable systems such as quadtrees and 
Delaunay refinement. We leave such developments to future research. 
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